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Abstract: 

We investigate the interface coupling between the 2D sine-Gordon equation and 
the 2D wave equation in the context of a Josephson window junction using a 
finite volume numerical method and soliton perturbation theory. The geometry 
of the domain as well as the electrical coupling parameters are considered. 
When the linear region is located at each end of the nonlinear domain, we 
derive an effective ID model, and using soliton perturbation theory, compute 
the fixed points that can trap either a kink or antikink at an interface between 
two sine-Gordon media. This approximate analysis is validated by comparing 
with the solution of the partial differential equation and describes kink motion 
in the ID window junction. Using this, we analyze steady state kink motion 
and derive values for the average speed in the ID and 2D systems. Finally we 
show how geometry and the coupling parameters can destabilize kink motion. 
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1 Introduction 



Nonlinear wave equations in one dimension that are close to integrable have 
been studied extensively for the past twenty years. The shift from exact in- 
tegrability can be due to the presence of extra terms in the equation or to 
an additional evolution equation. In the former case perturbation methods 
give good results. This kind of approach is far more difficult when two such 
equations are coupled. In that case one needs to resort to extensive numeri- 
cal studies in order to gain insight and develop approximate models. This is 
also true when one uses realistic boundary conditions, such as an impedance 
condition. 

Here we consider a 2D sine-Gordon equation coupled to a 2D wave equation via 
interface conditions. Such a model was originally derived for a so-called window 
Josephson junction [^,0 between two superconductors where the thickness of 
the oxide layer is small in the junction area and larger in the surrounding 
region. The electrodynamics of such a system is described by the sine-Gordon 
equation in the window area [Q where tunneling of electron pairs is possible 
and by the wave equation in the surrounding passive region. This model is 
much more general and could describe the motion of dislocations in a stressed 
inhomogeneous thin plate in this case the on site potential would be non 
uniform. Another example is the motion of domain walls in a one dimensional 
inhomogeneous ferromagnet One could also generalize the 2D rotator array 
of [0. From a fundamental point of view, the window junction allows to study 
in detail the interactions between a linear and a nonlinear system. The relation 
between the two subsystems can be changed by using a different geometry or 
topology or by modifying the coupling parameters. In the case of the Josephson 
junction, the former are the junction inductance and capacitance per unit area. 

Let us now recall what has been done to solve this problem. For static solutions 
the longitudinal extension of the passive region has been shown not to play a 
big role, only the lateral extension changes the solution 0. For small exten- 
sions w', the length of the problem is rescaled from Aj into X^s = \/l + 2^7xf 
1^ where Lj (resp. Lj) is the inductance in the junction (resp. passive region). 
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This rescaling explains very well the static solutions in the presence of a mag- 
netic field IP], it is a local effect which can also be seen for annular geometries 
For very large passive regions, the kink width is proportional to the inverse 
of the width of the junction, it can be very big and cause the destruction of 
the kink due to the longitudinal junction boundaries|I|. This could explain 
the absence of fluxon motion in window junctions with large passive regions. 

Lee and co-workers lTy,|TT]] investigated an infinite linear superconducting strip 
line where no Josephson current is present. They obtained the dispersion re- 
lation indicating that waves can propagate in the x direction and have a 
transverse structure in y. This study done in the linear case does not give 
any indications on the effect of the linear passive region on kinks which are 
strongly nonlinear. Another issue is that this analysis is based on eigenmodes 
and cannot be extended to the case when there is an external current applied 
to the device. 



In |]T^ experiments have been conducted on Fiske steps in window junctions 
with two different geometries a lateral passive region and longitudinal passive 
region. Lumped elements were assumed and simple models were derived from 
which the velocity was obtained. A passive region placed at each end of the 
junction acts as a lumped capacitance and generates radiation. On the con- 
trary a passive region placed along the junction acts as a transmission line 
in parallel to the junction. It gives a rescaling of the Josephson penetration 
length and a larger fluxon rest mass in agreement with Experimental re- 
sults for different extensions of the passive region together with preliminary 



numerical results have been reported in WM. These showed that fluxon motion 
ceases to be stable when the ratio of the widths of the passive region and the 
junction becomes larger than 3. The influence of the electric parameters was 
not studied. 

In a recent work one of the authors studied the case of a window junction 
with a homogeneous lateral passive region with periodic boundary conditions 
in the longitudinal direction. The motion of kinks in such a device occurs for 
velocities which can be calculated from the parameters of the device. In this 
study we will show that the presence of the passive region along the direction 
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of propagation will impose the velocity of the kink. 

We have studied the problem in a systematic way by comparing the solution of 
an effective ID problem with the 2D solution. We varied both the geometrical 
parameters and the electric parameters. In a first step we considered the prop- 
agation of a kink across an interface where the coefficient of the kink term and 
the electric properties vary abruptly. The operator describing this situation is 
discontinuous and one needs to use a finite volume approximation in order to 
have an accurate numerical solution. This method based on integrating the 
operator on reference volumes enables to satisfy exactly the jump conditions 
at the interface. Using this well suited numerical method together with soli- 
ton perturbation techniques, we have derived simple models explaining the 
dynamics of kinks in a window junction. 

The paper is organized as follows, in section 2 we introduce the 2D partial dif- 
ferential equation model. In section 3 we obtain an effective ID model for the 
case where the passive region exists only at each end of the junction. In section 
4 we study kink motion across an interface, we derive the perturbation equa- 
tions, analyze their fixed points and validate this approach by comparing the 
solution to the one given by a numerical integration of the partial differential 
equation. Section 5 discusses kink motion in a ID and 2D window junction. In 
section 6 we present kink motion instabilities in the system, consider limiting 
cases and give our concluding remarks. 



2 The model 



The electrodynamics of a window junction can be described by two 2D arrays 
of inductances coupled together through an RSJ element containing a capac- 



itor, the nonlinear Josephson element and a resistor] 15, 16 1. This approach is 
equivalent to a discretisation of Maxwell's equations together with Joseph- 
son's constitutive equation. We then assume perfect symmetry and go to the 
continuum limit to obtain the following partial differential equation for the 
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evolution of the phase difference in a domain Q 

C(f)u - div (;^V0) + e{x,y){sm{(f)) + a(f)t) = 0, (1) 

where C and L are respectively the normalized capacitance per unit surface 
and inductance and e{x, y) is the indicator function of the junction domain Qj 
(ie e = 1 in Qj and outside). The unit of space is the Josephson penetration 
depth and the unit of time the plasma frequency in the junction [Q]. 

Defining Cj and Lj to be respectively the normalized capacitance per unit 
surface and inductance in the passive region, Equation (|lD can be rewritten 
as the system 

^-A0 + sin0 + «^ = O in fij , (2) 



d'^^ 1 



together with the interface conditions for the phase and its normal gradient 
the surface current on the junction boundary dVLj 

^ = and — — = — (4) 
L/ on on 

where n is for example the exterior normal. The boundary condition on the 
boundary of the passive region represents the input of an external current or 
magnetic field 

Note that the jump condition (2nd relation in @j) can be obtained by inte- 
grating (|I]) on a small surface overlapping the junction domain Vtj. 

In the rest of the paper we have assumed a rectangular window of length / = 10 
and width w = 1 embedded in a rectangular passive region of extension w' as 
shown in Figure 1. We will not consider the influence of an external magnetic 
field and will assume the external current feed to be of overlap type so that 
the boundary conditions (13) become 
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1 dijj I , ,w 

L.Ty = ^W^) ^ = ±<2 

— — =0 for X = ±( — V w ). 
ox 2 



We assumed throughout the study a small damping a = 0.01 which is typical 
of underdamped Josephson junctions. 



3 The ID effective model 



In this section we introduce a simplified model where the passive region exists 
only on the longitudinal sides of the junction as in the bottom right panel of 
Fig. 1. Then the functions L and e depend only on the variable x. To simplify 
the problem we assume a uniform boundary condition ^ = T 2(i+2w') 
y = ±{^+w'). 

We will show that this system is well described in the limit of a small current 
/ by a one dimensional sine-Gordon equation. 



For that we write the solution of as <^ = + (pn where 

, ly' 

2w{l + 2w') 



(7) 



satisfies the y boundary conditions. These also imply that the residual tpn is 
even in y and can be expanded in a cosine Fourier series 



')R = J2^n{x,t) cos 



2mTy 



n=0 



w 



Inserting (|^ and (H) into (|I]) we obtain 



C[Aott + Ai^^ cos 



2iTy 



w 



1 / . . 277-0 , 

7(^0. + A. COS^ + ...^ 
L w 



1 
L 



ly 



w{l + 2w') 



2ti . 2Txy 
Ai sm h ... 

w w 



+e{x)[sm{(j)i+Ao+Ai cos —+...)+a{Ao^+Ai^ cos ^+...)] = 



w 



w 
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We get the evolution of Aq by integrating the equation above in y and dividing 
by w, 



CAo,-(-AoJ 4 



+e{x) 



1 2 

— / sin(07 + + cos \- ...)dy + aAo 

w J w 



We write 

^TTy ^Tiy 

sin(0/ + Aq + Ai cos h ...) ~ sin(0/ + Aq) + Ai cos cos(07 + Aq) 

w w 

because 1^41 1 << \Aq\. First we consider the integral of the first term sin(0/ + 

Aq) 

2 2 2 

sin(07 + AQ)dy = cos(Ao) dysn\{4)i) + sin(>lo) ^ cos(0/) 

— w — ttJ — w 

2 2 2 

The first integral on the right hand side can be written 

2 1 

I dysvaidr) — w dCsin( -C'^) 

J ^ ' J ^ ^ 8 / + 2w' 

-w Q ^ ' 

2 



SO that if — g(-;^2w') 1 the sine can be linearized and the cosine taken equal 



to 1 leading to 



Iw 

sin((/)7 + AQ)dy = - 24.(1 + 2w') ^^^^^^^ + sin(Ao) . 

— w 
2 



The terms in Ai are 



dy cos{4>i) cos sin(/lo) J dy sin(0/) cos 



27ry 
w 



Both terms can be treated following the same approximations as above so that 



dy cos{(f)i) cos 



27rv w 



W TT 



dy sin(0/) cos 



2'Ky Iw"^ 
w ~ AiT^{l + 2w') 
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Combining all contributions, we obtain for the evolution of Aq 

Iw w Iw^ 

- cos(Ao) + cos(Ao) - .3..^^^ sin(Ao)) + «AoJ = 0. 



24(/ + 2w') ' V ' 47r3(/ + 2w;' 

For the values of the geometric parameters that we have taken / = 10,tf = 
1, \ w'\ < 10 and if we assume I < 1 which is the case for a zero field step then 
the terms 24{i+2w') ' f ' 4t7^i+2w') neglected so that we are left with the 

following sine-Gordon equation for Aq (dropping the 0) 

CAu - Qa,) + 7(x) + e(x) [sin(A) + a A] = 0, (10) 

where 7(x) = i(^^)^l^i^2w') - '^^^^ approach is validated for a homogeneous 2D 
sine-Gordon equation by the numerical evidence provided by Eilbeck and al 



23 that when {I/8C ^ 1), the phase is almost uniform in the y direction. 



4 Motion of a kink across a ID interface 



Before considering the motion of a kink in a window junction, we will study the 
simpler problem of an interface between two media with different values of the 
parameters {L{x), C{x), e}. The propagation of a soliton across an interface 
has been studied by many authors lH], |2|, A pioneering work was done 



by Aceves et al for the motion of a nonlinear Schroedinger soliton between 
two media of different Kerr indices. It was shown that the soliton perturbation 
equations give a qualitatively correct description even when the nonlinearity 
of the second medium is very small. We will follow the same route and derive 
adiabatic equations for the soliton parameters. These are in principle only 
valid for small perturbations but some ideas can be obtained by taking the 
limit e ^ 0. In this way one can understand qualitatively the motion of a kink 
from a junction to a passive region. 

Following [|l^, we introduce the change of variable 



dz = JL{x)dx , (11) 
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which makes the perturbation theory regular in the sense that the solution 
has the same form on the left and right of the interface. In principle the terms 
CAtt and e sin A make the perturbation theory singular, however we are mostly 
interested in the fixed points that exist due to the interfaces, The first term 
will not influence these, nor their stability. We expect this description to hold 
if the wave remains close to the interface. 

Then we write equation (p!OD in standard form 

Att - A,, + sm{A) = - 7(s) - aAt + (1 - e{x)) sm{A) + (1 - C{x))Au 

- [InL(x)] A, . (12) 

To simplify the calculations we assumed a uniform damping a. To model the 



experimental situation of [Tj] for example, it is natural to assume the following 




dependencies for the parameters L, C and e 



lifx<0 lifx<0 
L{x) = \ , C(x) = < and e{x) 

if a; > if x > 



With this type of distribution of the inductance L and using (|Tl|) we obtain 
the differential relation [InL(x)] = {\\i{Li))5{z) , where 5 is the Dirac delta 
function. 

Now we rearrange the right hand side of Eq. (|12D as the perturbation term 

-oyz) 

(13) 

where /iq = cq - 1, 7/ = ^^^(/_^2^') > 7 J = ^ ^^^^^ Heaviside 

function. 

Note that the functions C and e are piece-wise constants and therefore there is 
no different re-scaling on each side of the interface as we go from the variable 



ef = -Ij - aAt + [-(7, - 7^) - ^0 sin A + {1- Ci)Att] H{z) - ^^5{z)A 



X to z. Having set up the problem ([T^) , we assume as usual that the soliton 



parameters are slowly modulated and use perturbation theory to derive 
their evolution. 
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We choose the kink ansatz 



Aa{z,t) = A tan ^expa^ 



z-Z{ty 

7r^, 



where Z{t) = J v{t')dt' + Zo{t) and a = 1 (resp. —1) for a kink (resp. an- 



tikink). 

Then {Z{t),v{t)) are solutions of the differential system 



dt 2 I 



X sech 



(1 - CiyvT^ 



- V' 

Z 



1 + tanh 



-(7/-7j)Vr^2G 



(14) 



2^l 



dv _ -ncy^^i + 7j) (1 - v ) 
~dt ~ 8 

X sech^ 



av{\ — V ) — 



2, 



z \ 



- ln(L, 



(15) 



where 2G = 1.831931188.. [0. If we assume e = 1, a = / = and = 1 



these equations are exactly the ones derived in |T9 



In the absence of the current and the damping, the equations (|n[|T^) derive 
from the Hamiltonian, 



H 



+CX3 



C{z) (dv\\l(dv\\ , ■ 



dz . (16) 



If we assume the fiuxon to have velocity V[ for z -C 0, its energy is Hi{vi) = 
(i_,|)i/2 • If it crosses into the passive medium 2; 0, it will have a velocity Vr 
and energy Hyivr) 



Hr{Vr) 



Avl{Ci-eo)+A{l + ep) 



(17) 



and one can obtain the condition on the velocity vi so that the left coming 
soliton crosses into the right hand medium as Hi[vi) > Hr{vr = 0) which 
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implies vf > 1 — ■ One can then compute Vr by identifying Hi and Hr- 

If vi < ^1 — (i^~)2 the fluxon is reflected with velocity —vi. 

4-1 Existence and stability of fixed points 

From the above equations, we obtain two fixed points symmetrically placed 
with respect to the interface z = 



Z± = atanh rj^ = atanh ± [ 1 ' 



Co - 1 - log(L/) J 
where 7 = 

The expression indicates that at least two ingredients are necessary for soliton 
trapping at a finite distance from the interface, a DC current 7, a jump in 
the coefficient of the sine nonlinearity e (the critical current density in the 
Josephson model) and an inductance jump Lj. Trapping has been predicted 
with the first two features in the works pO|j2T[| . Notice also that the existence 
and position of the fixed points is completely independent of the damping a 
and capacity C/ in the passive region. 

Let us consider the conditions for existence of these fixed points. For that we 
separate the cases of a kink or an antikink. For a kink o" = 1, the fixed points 
exist if fiQ — log{Lj) — 717 > or L/ < e^°~'^'^ while for the antikink a = —1 
they exist if Lj > 6^°+'^'^. 

To investigate the stability of these fixed points we compute the Jacobian of 
the evolution equations ( p!^ - |T5| ) and estimate its eigenvalues A at each of the 
fixed points. The characteristic equation is 

X +Xa -— [1 - —(1 + 7]±) - -Z±-f(x — ^ = 0. 

2 2 2 jjiQ — log L/ 

Let us fix the type of wave to be an antikink. Then we find that the discrimi- 
nant for the fixed point Z+ is negative for any value of the current and tends 
to zero as the current tends to zero. This fixed point is then stable because 



11 



Fig. 1. Table 1: existence and stability of the fixed points Zj^ 



Li < e^o-^-'^T 



g6o-l+7r7 < 



antikink [a — —1) 



no fixed points 



no fixed points 



2 fixed points Z± 
is stable 
Z_ is unstable 



kink (a = 1) 



2 fixed points Z-t 



Z^ is unstable 



Z_ is stable 



no fixed points 



no fixed points 



the eigenvalues are complex conjugate and their real part —a/2 is negative. 
On the other hand the fixed point Z_ is unstable because the discriminant is 
always positive and the eigenvalues are of opposite signs. 

If we had considered as for the 2D problem that dissipation is absent on 
the right hand medium, we would have obtained half the damping term in 
the equation for dv/dt and an additional term — v'^v^alog{2) in the 

equation for dZ/dt. The first term would affect the stabihty, it is equivalent 
to halving the overall damping term. 

The existence and stability for both kink and antikink waves is summarized 
in Table 1. 

This study that such an interface enables to trap a certain type of wave. 
For large L/, antikinks are trapped at a fixed point X+ — Z^/^Lj) in the 
passive region while kinks remain free. On the contrary when Lj is small, 
kinks are trapped at the fixed point X_ — Z_. If the interface is now such 
that the passive region is for 2; < and the junction for 2; > 0, then the 
second column and the last column of the table should be permutated and 
the stability properties exchanged. This provides a qualitative understanding 
of the window device which contains two such interfaces. 
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4-2 Validation of the perturbation theory 



In order to test these predictions we have compared the solution of the par- 
tial differential equation (TU) with the solution of the perturbation equations 



(p!^ - |15D . Note that the former is a generalized operator with discontinuities 
in the coefficients so that the solution is continuous but exhibits jumps of 
its derivative at the junction interface, To obtain it we introduce the finite 



volume method used for the treatment of hyperbolic equations The idea 
is to integrate the partial differential equation over intervals around each dis- 
cretisation point where the solution is assumed constant. In this way the jump 
conditions are exactly respected. The details of the implementation for both 
the ID and 2D cases are given respectively in appendices A and B. 

Fig. 2 shows the phase plane in the original position- velocity coordinates {X, u) 
for an antikink wave for two different interfaces with eg = 0.6,7 = 0-02 and 
L/ = 1 (top panel) and L/ = 1.5 (bottom panel). The solution of the perturba- 
tion equations is given in full line while the solution of the partial differential 
equation is given by the crosses. A least square procedure has been used to 
estimate the position and velocity of the wave. For both cases we obtain a 
good overall agreement between the solution of the partial differential equa- 
tion and the perturbation approach, showing that the fixed point traps the 
antikink even for a large initial velocity [u = 0.5). For the bottom panel where 
Lj = 1.5 notice the jump in the phase gradient and the fact that the fixed 
point is close to the interface as expected. 

We now consider that the nonlinear term on the right hand side is very small. 
Fig. 3 shows X{t) and u{t) for an antikink in the case Lj = 1 and eo = 0.1. 
For such a strong perturbation, we obtain a qualitative agreement for both the 
position and velocity of the wave between the solutions of (|l^) and ([T^ITsD . 



Both the numerical solution and the perturbation method show the existence 
of a stable fixed point located in the region e(x) = cq < 1. We therefore expect 



([T^ - |15D to provide qualitative results even for cq ^ i.e. the nonlinear-linear 
interface. It is then possible to use these arguments to describe the motion of 
a kink in a window junction where the passive region exists on both sides of 
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the Josephson junction. 



5 Motion of a kink in a window junction 



5. 1 The perturbation theory 



As shown in the perturbation equations ( p!4^T5|) the motion of a kink is due 
to the injection of direct current in the device, described by the 7 term. In 
a pure Josephson junction, this gives rise to the well-known Zero (magnetic) 
Field Steps in the current voltage IV characteristics. The presence of passive 
regions where no sine term is present and for which there is an inductance 
jump will affect the kink motion. We will concentrate in this study on the 
features of this motion. The study of the IV characteristics and its associated 
zero field step will be presented in |T7|. 



The study conducted above for the case of a single interface provides a way 
to give a simple description of the fiuxon motion in a one dimensional window 
junction (bottom left panel in Fig. 1). In that case the current density is 
lidw = -^(1 + '^1 Li)/{2{1 + 2w')) and the position of the fixed points is 

X± = ^atanh ± f 1 + (T ^ V ' ^^^^ 

Li \ eo - 1 - log(L/) / 

This quantity depends weakly on both the current / and the inductance in 
the passive region L/. For the parameters used here eo = 0, / = 10, the fixed 
point is inside the domain as soon as w' > 2. 

Fig. 4 shows the schematic motion of a wave in such a device, in the limiting 
cases Li » 1 (top panel) and Li << 1 (bottom panel). The antikink wave 
propagates towards the right and gets refiected as a kink at the right boundary 
of the device. When L/ >> 1, the antikink finds two fixed points on each side 
of the interface, identified by large or small circles depending whether the fixed 
point is stable or unstable. The antikink is slowed down and can be trapped 
at the fixed point located to the left of the left interface. The kink can become 
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trapped at the fixed point of tlie riglit interface. The situation is reversed for 
Li « 1. We see that these points will affect the motion and could hinder 
wave motion. 



5.2 Numerical study 



To confirm the role of the fixed points at the interfaces we carried out a sys- 
tematic study of the kink motion in both a ID and a 2D configuration. The 
numerical procedure in both cases is to discretize the spatial operator using 
the finite volume approximation which allows to preserve the interface con- 
ditions (see Appendix A and B). The temporal part is then advanced using 
the Dormand and Prince ordinary differential equation solver D0PRI5 imple- 
mented by Hairer and Norsett [^. The initial condition is a static kink for 
the ID junction and the static solution [0 in the case of a 2D window junc- 
tion. The wave is then accelerated by the injection of direct current through 
the boundary conditions (P). It will reach its limiting velocity which can be 
obtained by simple arguments. 

Figure 5 presents the motion of a kink at its limiting velocity in a ID window 
junction. The top and middle panels show respectively the position X{t) and 
velocity u{t) vs. time while the bottom panel shows the phase plane (X, -u). 
The electric parameters are L/ = C/ = 2 so that the velocity of the linear 
waves in the passive region is vj = l/y/LjCj = 0.5. One sees that the kink 
velocity is close to 1 in the junction region while it is about 0.5 in the passive 
region, so the kink adapts its speed to the region it travels in. At this point 
it is interesting to remark that the kink continues to exist even in the regions 
where there is no sine term. The whole motion and the voltage observed are 
then given by the expression 

A0 27r 

The motion of a kink in a 2D window junction is different as can be seen from 
Figure 6 which shows a contour plot of the phase ip{x,y,t) for 4 successive 
values of time. The parameters are the same as in Figure 5. In this case 
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no significant change is observed in the speed of the wave in the junction and 
passive regions. The voltage observed indicates that the velocity is everywhere 



equal to the one in the passive region vj = l/y^LjCj. There is a simple 
explanation for this, the soliton velocity is a free parameter for the sine Gordon 
equation. In the linear region waves can only travel at velocity vj. The dressed 
kink which carries a significant part of its energy in the passive region adapts 
its speed to vj. The voltage in this case is 

= (21) 
2w' + l ^ ^ 

For a large extension w' of the passive region or values of the electric parame- 
ters non equal to 1, the kink motion can break down. In the next section, we 
give a few examples of this phenomenon. 



6 Instabilities and limiting cases 



In our procedure, we take as initial condition a static kink and accelerate 
it by injecting current via the boundary condition (^. We present now a 
few situations where this procedure led to the break down of the kink and 
another type of dynamical behavior occurred. The situation is simpler for the 
ID window junction and we will present this first. 

In the ID case, the kink (antikink) can be trapped at the fixed point on the 
right (left) end interface. This occurs because as the width of the passive region 
w' is increased, the current density is decreased so that there is less driving 
force to overcome the potential barrier created by the stable fixed point. For 
the junction geometry that we chose I = 10,w = 1 kink motion breakdown 
occurs for w' ~ 10. 

The electrical parameters act differently. To see this consider the equations of 
the problem 

_|__| + „„^ + , + „_| = N<-, (22) 
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together with the interface condition at |a;| = | 

, , 1 dip d(j) 
'^ = ^ ^U^ = d-x ' 

and homogeneous boundary conditions ^\x=±{i/2+w') = 0. 

When Ci or L/ are large we can consider their inverse as a small parameter 
and use perturbation theory to get some estimates. Let us first consider the 
case C/ ^ 1. Then from equation (^) we get ^ = so that ^ = is a 
constant. We then write ip = Vt + f{x) so that the ( P5| ) becomes 



+ 7^/ = 



which yields f- = •yLj^x - 1/2 - w') for 1/2 < x < 1/2 + w' and ip = •yLj^x - 



dx 



1/2 — w') /2 + Vt + Cte. This fixes the boundary conditions for the first 
equation. When Li is large ^\x=±i/2 — > so that kink motion is possible. On 
the contrary, kink motion becomes impossible when Lj << 1 because then 

0^\x=±l/2 OO. 

Let us now consider the two dimensional situation. For large passive regions 
the static kink solution becomes stretched and occupies the whole junction 
It is then difficult to accelerate and generally it breaks up after a few round 
trips. Figure 7 shows contour plots of the phase ip{x,y,t), numerical solution 
of (|1]) for t = 0, 5, 6 and 8 clockwise starting from the top left panel, for a 
large passive width w' = 7. One sees that the kink breaks up and leads to a 
radial phase distribution centered on the junction. The wave equation then 
dominates the dynamics and the average phase increases with time. 

Another type of instability is due to the electrical parameters Lj and Cj, 
especially the former because it acts on the interface condition in addition 
to the operator. For example consider Lj = 10~^ and w' = 2, the case of 
Figure 8 where the phase (p{x,y,t) is shown for t = 2, 10 and 20 from left to 
right. The top panels show the contours while the bottom ones show the three 
dimensional plots. The phase tends to obey a Laplace equation in the passive 
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region and there the kink leads to a uniform phase distribution. The junction 
is strongly coupled to the passive region because the gradient at the interface 
is very large. This causes the break-up of the initial kink and leads to a radial 
solution with small oscillations. 

A very large inductance leads to the creation of a boundary layer along the 
top and bottom edges of the device as shown in Figure 9 which presents two 
snapshots for successive times. In this case L/ = 10^. At the interface the 
gradient of if is close to zero, so that the phase is constant in the junction. 
Notice the strong gradients on the top and bottom sides of the passive region. 

For w' = 2, we have found in the [Lj, Cj) plane, the region of instability of the 
kink as shown in Figure 10. We plotted in the (L/, C/) plane, the hyperbolas 
corresponding to three values of velocity in the passive region, vj = 0.5, 1 and 
2. The signs (*) (x) and (+) indicate where we found a stable kink motion. 
We have isolated by the solid curve the stability region in the [Lj , C/) plane. 
As expected the instability is due essentially to the inductance Lj so the kink 
is stable in the domain 0.1 < Lj < 2. In this interval we always observed 
kink motion independently of the value of the capacity Cj. For example for a 
capacity Cj = 10^ and Lj = 1, we obtain a stable kink motion. The voltage 
observed Vnum2-D = 4.510 x 10~^ is in good agreement with the one given by 
(|21|), V2D = 4.5 X 10~'^. Then we conclude that the region of stability of the 
kink in the window problem is 

^stability = {0.1 <Lj< 2}. (24) 

At this time we do not have an explanation for this treshold value Lj = 2. 

From a general point of view this study has shown the importance of the per- 
turbation approach to gain insight into a complicated wave problem. Another 
determining factor was the use of the finite volume discretisation to solve the 
inhomogeneous partial differential equation. 
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APPENDIX A: The finite volume formulation in ID 

We introduce new functions for the ID window equation (p!0| ) 

u = At and v = , 



which becomes the system 
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where 





u 




u = 




, A{x) = 




V 









1 

L{x)C(x) 

1 



and 



B{U,A) 



tlx) T 

^ ' (aw + sin A) - ' 



C{x) 



C[x) dx \C{x) J L{x) 



u 



The matrix A has two real eigenvalues Ai,2 = ^ y l{x)c(x) ^'^^^ system 
p5| ) is strictly hyperbolic with a source term. Due to the discontinuities of 
the functions C{x),L{x) and e(x) at the interfaces, the numerical integration 
of the Eq. ( ^5|) using finite differences cannot be done. We therefore turn to 
the finite volume method |]2B|, whose principle is first to consider a partition 
of the passive and junction domain in cells 



x\Xk Y - ^ ^ 2 



where Xk are collocation points defined by 



Xk = {k — in passive region , 
Xk = {k — l)hj in junction region . 



Here hi and hj are the space mesh-length, defined by 



w 



hi 



I 



where rij and rij respectively are the number of discretization in the passive 
and junction domain. The total number points of the whole domain is 2ni + nj. 
Notice that the points Xk such k = rii and k = rii + rij correspond to the real 
interface points x = —1/2 and x = 1/2. 

In a second step we assume the solution of Eq. (pSj) to be constant in each cell 
Vk''- We integrate Eq. (p5|) over each cell 



9U I djAU) 
dt dx 



Bdx, 



(26) 
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and obtain 



where Uk = U{x) for all x e V^'^ . 

We therefore obtain the following system of ordinary differential equations In 
the passive region 



d 
di 







Vk 





{vh{Xk + I) - Vh{Xk - f 
u-h {^k+^^ -Uh i^k-^^ 



which in terms of A yields 

dt 

duk _ 
dt 





7 




Ci 








Ak+1 - 2Ak + Ak-i 
Cj LjCjhf 



(27) 



(28) 



In the junction domain we obtain 



d_ 
di 







J_ 


Uk 




hj 


Vk 







sin Lpk + aUk + 7 




, (29) 



which in terms of A implies 
dAk 



dt 



Uk 



duk . . , Ak+i - 2Ak + Ak-i 
— = -aUk - sm - 7 H 



(30) 
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At the interface we get 





Ci 




1 






Uk 









Uh (xk + ^) 


^ h, 
2 


1 


2 


1 


J 


d 


Vk 




1 




Vh {xk + ^) 



1 



1 



Vh\Xk- ^ 



which in terms of A yields 

dAk 



dt 



= Uk 



duk 



A 



k+l 



hi 



\ {hi + hj) + -f{auk + sin Ak) 



Ak Ak - A, 



k-l 



(31) 



Lihi 



hj , . . X 1 hi + hj 

y [auk + sm Ak) - 7 I — ^ 



APPENDIX B: The finite volume formulation in 2D 



Here we follow a similar procedure as in the ID case. The discretization of 
the domain in the x direction is the same as done in the appendix A and in y 
direction we introduce the fellowing partition 



y\yk-^ <y<yk + ^ 



where i/k are the collocation points defined by 



yk — ~ 7:)K passive region , 



yk — {k — i)h', in junction region . 



Here /i^ and h'j are the space mesh-length along y direction, defined by 



' K + y^ n;.-i' 
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where n[ and n'^ respectively are the number of discretization points in the 
passive and junction domains. The total number points in the whole domain 
is {2ni + rij) x (2n- + n'^). 

For convenience we choose the same step size in the whole domain (the number 
of points ni,nj,n^,n'j are chosen such hi = hj = h[ = h'j = h). Then the 
discretization of the domain can be written as 

= [J VLk^k' , 

k,k' 

where 

^k,k' = l^{x,y)\xk-^ < X < Xk + ^ and Vk' - ^ < y < Vk' + , 

are squares of side lenght h centered at the point {xk,yk')- 

We integrate equation (0) on each cell flk,k' and use the finite volume approx- 
imation to obtain 

i^k,k' J C dxdy- J -^V0-ndcr+(aVfc,fc'+sin(0fc,fc/)) J e dxdy = , (33) 



where ■0 = 0, n is the unit exterior normal to dflk,k' and (pk^k' = 4>{x,y), 
'^k,k' = 4't{x,y) are assumed constant for all {x,y) G ^lk,k'- We then obtain for 
a cell inside the junction 

4>k,k' = '^k,k' 

'>Pk,k' = ^[4>k+i,k' + 4>k,k'+i + 4>k-i,k' + 4>k,k'-i - 40fe,fc'] + sin (j)k,k' + Oi'ipk,k'- 

(34) 

In the system defined above M, M' respectively are the total number of points 
in the x and y directions. For k = 1,M and k' = 1, M', the solution {(pk,k' , 'ipk,k') 
is determined using the boundary condition (|^). 

For a cell in the passive region we obtain 



4'k,k' — '^k,k' 

1pk,k' = h'^LiCi i^k+l,k' + (t>k,k'+l + (pk,k'-l — '^<t>k,k')- 
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Now consider cells that are overlapping the interface. We will dislpay two cases 
to show how the method works. First consider the cell overlapping the 
bottom left corner of the junction boundary (x^ = — //2,?/^ = —w/2). 
There ( |5B| ) yields 



-M'|(i7 + i) + <^M'+iKrF + i) + 



Second consider a cell that is overlapping on the bottom boundary of 
the junction {—1/2 < Xk < i/2,y'i^ = —w/2). In that case we obtain 



)k,k' 



k,k' 



k,k' 



>k-l,k' 



fe^(Cj+i) [0fc+i,fc'2(rF + 1) + 4>k,k'+i + 

- 0fc,fc'(2 + ;^)] + Y:f^(sin0fe,fc' + ai)k,k' 



With this formulation we have checked the conservation of energy in the ab- 
sence of current or damping a = / = 



E 



f(^) +^|V^P + (l-cos(0))]dxdy. 



For a number of points in x and y, M = 100, we obtain that the energy is 
conserved up to 10~^ in relative error for t < 1. 
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FIGURE CAPTIONS 



1 Top panel: a view of a window Josephson junction. The bottom left 
panel shows a schematic top view. For the system shown on the right the 
linear region exists only on the left and right sides of the junction. 

2 Phase plane {X, u) for an antikink wave propagating across an interface 
at a; = 0. The solid line presents the solution of the adiabatic equations 
([I^plSD and the losanges are the numerical solution of the partial differ- 
ential equation (|T0|). The parameters are eo = 0.6,7 = 0.02, L/ = 1 (top 
panel) and L/ = 1.5 (bottom panel). 

3 Interface between a nonlinear and linear medium, the top panel shows a 
plot of the antikink position X{t) vs. time for the solution of the adiabatic 
equations (|l^|l5]) in full line. The crosses indicate the solution of the 
partial differential equation (|10|). The bottom panel presents the velocity 
v{t). The parameters are eo = 0.1,7 = 0.01, L/ = 1. 

4 Simplified description of the fluxon motion in a one dimensional win- 
dow junction for L/ ^ 1 (top panel) and Li <^ 1 (bottom panel). The 
stable fixed points are shown by the large circles and the unstable ones 
are the small circles. 

5 Motion of a fluxon in a one dimensional window junction from the 
numerical solution of ([T0|). The top and middle panels show respectively 
the position X{t) and velocity u{t) vs. time while the bottom panel shows 
the phase plane (X, m). The parameters are eo = 0,w' = 2,Li = Cj = 
2,1 = 0.1.. 

6 Contour plot of the phase ip{x,y,t) for t = 0,400,600,800 of the nu- 
merical solution of the two dimensional window junction equations (|I]). 
The parameters are the same as in Figure 5. 

7 Contour plots of the phase ip{x,y,t), numerical solution of (|l|) for t = 
0, 5, 6 and 8 clockwise starting from the top left panel. The parameters 
are w' = 7, Lj = Ci = 1, 1 = 0.3. 

8 Plots of the phase ip{x, y, t), numerical solution of (P for t = 2, 10 and 
20 from left to right. The top panels show the contours while the bot- 
tom show the three dimensional plots. The parameters are w' = 2, Lj = 
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0.0001, = 1,1 = 0.3. 

9 Contour plots of the phase ip{x,y,t), numerical solution of (|l|) for t = 
50100 (top) and 51000 (bottom) from left to right. The parameters are 
the same as in Fig. 8 except Lj = 10000. 

10 Parameter plane {Lj, Cj) showing the regions of existence of zero field 
steps (ZFS) corresponding to the shuttling motion of a fluxon. The ve- 
locities vj are indicated in the top right corner of the picture. 



27 




Fig.l , Benabdallah et al 




Fig. 2 , Benabdallah et al 




Fig. 3 , Benabdallah et al 



Large inductance 
antikink 



wave 
equation 


Sine-Gordon o 


wave 



equation 


kink - — ^/ 


wave 
equation 


Sine-Gordon 


wave 
equation 




small inductance 






antikink V — ► 




wave 
equation 


Sine-Gordon 


wave 
equation 


kink - — ^ 


wave 
equation 


Sine-Gordon o 


wave 

o 

equation 



4 , Benabdallah et al 
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